%% Navier-Stokes equation solver
%
%
%

%% Problem settings

clear all
clc

global NIt
global NJt
global NIe
global NJe

% Domain dimensions

Lxt  = 12;
Lxe  = 4;       % Longitude before expansion
Lyt  = 2;
Lye  = 1;

% Discretization, pressure finite volumes

% I denotes x direction
% J denotes y direction

elemsPerLenghtUnit = 5;

NIt = Lxt*elemsPerLenghtUnit;
NIe = Lxe*elemsPerLenghtUnit;
NJt = Lyt*elemsPerLenghtUnit;
NJe = Lye*elemsPerLenghtUnit;

% Nodal spacing

h   = Lxt/NIt;

% Constitutive parameters

nu  = 0.01;
rho = 1;

% Boundary conditions


% Time domain and discretization
dt = 0.1;
T0 = 0;
T1 = 1.;
nSkip = 1;             % skips saving solution
nTimeSteps = (T1-T0)/dt;
nSaveSteps =  (T1-T0)/(dt*nSkip);
nSave      = 0;

dtLim1 = 0.25 * h^2 / nu
dtLim2 = 2*nu/1

if (dt>dtLim1 || dt>dtLim2 )
fprintf ('stability warning')
end



%% Pre-processing
%  Memory allocation

dimP = NIt*NJt-NIe*(NJt-NJe);
dimU = dimP + NJt;
dimV = dimP + NIt;

% Problem variables
P = zeros(dimP,1);
U = zeros(dimU,1);
V = zeros(dimV,1);

% Auxiliars: midstep variables & Pk for iterative pressure scheme computation
Pk = zeros(dimP,1);
Us = zeros(dimU,1);         % U star
Vs = zeros(dimV,1);
Ua = zeros(dimU,1);         % Previous step solution
Va = zeros(dimV,1);

% Saving solution over time
PT = zeros(dimP,nSaveSteps);
UT = zeros(dimU,nSaveSteps);
VT = zeros(dimV,nSaveSteps);

% Auxiliar linear matrix & load vector for pressure implicit equation
M = zeros (dimP,dimP);
l = zeros (dimP,1);


% Dimensions

NIut = NIt+1;
NJut = NJt;
NIue = NIe;
NJue = NJe;

NIvt = NIt;
NJvt = NJt+1;
NIve = NIe;
NJve = NJe+1;

% Mapping spatial points

Xp = PspatialMap( h , Lxt , Lxe , Lye , Lyt , NIt , NJt , NJe , NIe , dimP);
Xu = UspatialMap( h , Lxt , Lxe , Lye , Lyt , NIut , NJut , NJue , NIue , dimU);
Xv = VspatialMap( h , Lxt , Lxe , Lye , Lyt , NIvt , NJvt , NJve , NIve , dimV);

assemblePressureMatrix;
% Vectorial implementation
assembleDMatrix;
Dp = sparse(D - (2*nu/dt)* eye(dimU+dimV));
dirichletConditionsDp;
Dpp = sparse(-D - (2*nu/dt)* eye(dimU+dimV));

% A_h( u^{n})
Au_vec = zeros(dimU+dimV,1);
% A_h( u^{n-1})
Aum1_vec = zeros(dimU+dimV,1);
%% Processing - Predictor step
%
%
%
%

% Temporal loop

time = T0;
nstep = 0;
while (time<T1)
    
    %if (nstep>1)
    Um1 = [U; V]; % vector [u v] at previous time
    %end
    
      
    
    time = time + dt;
    nstep = nstep+1;
    
    fprintf('Solving timestep %d of %d \n',nstep,nTimeSteps)
    
       
    % Velocity in X
    for i=1:NIut
        for j=1:NJut
            
            if (i<=NIue && j>NJue) continue; end
            
            U_cen   = U ( Umap (i  ,j  ) );  % center of finite volume
            
            if (j<=NJue)
            
            % Boundary conditions on inlet & outlet
            if (i==1) Us ( Umap (i,j) ) = getParabola ( Xp ( Umap (i,j),2 ) ) ; continue ; end    % Boundary condition inlet
            if (i==NIut) Us ( Umap (i,j) ) = Us ( Umap (i-1,j) ) ; continue ; end    % Boundary condition oulet - null normal gradient
            
            else
                
             % Boundary conditions on expansion & outlet
            if (i==NIue+1) Us ( Umap (i,j) ) = 0 ; continue ; end    % Boundary condition expansion
            if (i==NIut) Us ( Umap (i,j) ) = Us ( Umap (i-1,j) ) ; continue ; end    % Boundary condition oulet - null normal gradient
            
            end
            
            
            Ax = getAx( i,j,U,V,NJue,NIue ,NJut,h) ;
            Dx = getDx( i,j,U,V,NJue,NIue ,NJut,h) ;
            
            Us ( Umap (i,j) ) = U_cen + dt* (-Ax+nu*Dx);
            
            
            % vectorial ---------------------------------------------------            
            Au_vec(Umap(i,j)) = Ax;
            Aum1_vec(Umap(i,j)) = getAx( i,j,Um1(1:dimU),Um1(dimU+1:end),NJue,NIue ,NJut,h);
            
        end
    end
    
    
    % Velocity in Y
    for i=1:NIvt
        for j=1:NJvt
            
            if (i<=NIve && j>NJve) continue; end
            
            V_cen  = V ( Vmap (i  ,j  ) );  % center of finite volume
            
            if (j<=NJve)
            % Boundary conditions on top and bot boundaries
            if (j==1) Vs (Vmap (i,j)) = 0; continue ; end
            if (j==NJve && i<=NIve ) Vs (Vmap (i,j)) = 0; continue ; end
            else
             % Boundary conditions on top and bot boundaries
            if (j==NJvt ) Vs (Vmap (i,j)) = 0; continue ; end
            end
            
            Ay  = getAy( i,j,U,V,NJve,NIve ,NJvt,h);
            Dy  = getDy( i,j,U,V,NJve,NIve ,NJvt,h);
            
            Vs ( Vmap (i,j) ) = V_cen + dt* (-Ay+nu*Dy);
            
            % vectorial ---------------------------------------------------
            indexV_cen = Vmap(i,j);
            Au_vec(dimU+indexV_cen) = Ay;
            Aum1_vec(dimU+indexV_cen) = getAy( i,j,Um1(1:dimU),Um1(dimU+1:end),NJue,NIue ,NJut,h);
            
        end
    end
    
    
    % Mount vector
    U_vec = [U; V];   
    % RHS = (3/2)*(Dt) A_h( u^{n}) - (1/2)*(Dt) A_h( u^{n-1}) - D_h(u^n)
    % - (2\nu)/(Dt) u^n
    RHS = (3/2)*(dt)*Au_vec - (1/2)*(dt) * Aum1_vec - Dpp * U_vec;
    
    % RHS = BOUNDARY CONDITIONS 
    RHSboundaries;
    
    
    Ustar_vec = Dp \ RHS;
    
    Us = Ustar_vec(1:dimU);
    Vs = Ustar_vec(dimU+1:end);
    
    
    %% Processing - pressure step
    %
    
    %c2 = rho/(dt*h);
    c2 = 1/(dt*h);
    
    for i=1:NIt
        for j=1:NJt
            
            if (i<=NIe && j>NJe) continue; end
            
            U_right = Us ( Umap (i+1,j  ) );
            U_left  = Us ( Umap (i  ,j  ) );
            V_top   = Vs ( Vmap (i  ,j+1) );
            V_bot   = Vs ( Vmap (i  ,j  ) );
            
            if (i~=NIt)
                l ( Pmap (i,j) ) = U_right - U_left + V_top - V_bot;
            else
                l ( Pmap (i,j) ) = 0 ;
            end
            
        end
    end
    
    l = l/c2;
    %     if (nstep==1)
    %         M(1,:) = 0;
    %         M(1,1) = 1;
    %         l (1) = 0;
    %     else if (nstep==2)
    %             assemblePressureMatrix;
    %         end
    %     end
    
    P = M\l;
    
    
    
    %% Processing - projection step
    
    c1 = (dt/(h));
    
    % Velocity X
    
    for i=1:NIut
        for j=1:NJut
            if (i<=NIue && j>NJue) continue; end
            
            % Top boundaries
            if ( j==NJut ) U ( Umap (i,j) ) = Us ( Umap (i,j) ) ; continue ; end
            if ( (j==NJue) && (i<=NIue) ) U ( Umap (i,j) ) = Us ( Umap (i,j) ) ; continue ; end
            
            % Bot boundary
            if ( j==1 ) U ( Umap (i,j) ) = Us ( Umap (i,j) ) ; continue ; end
            
            % Right boundary
            if ( i==NIut ) U ( Umap (i,j) ) = U ( Umap (i-1,j) ) ; continue ; end
            
            % Left boundary
            if ( i==1 ) U ( Umap (i,j) ) = Us ( Umap (i,j) ) ; continue ; end
            if ( i==NIue+1 && j>NJue ) U ( Umap (i,j) ) = Us ( Umap (i,j) ) ; continue ; end
            
            
            
            U ( Umap (i,j) ) = Us ( Umap (i,j) ) - c1 * ( P ( Pmap(i,j) ) - P ( Pmap(i-1,j) ) ) ;
        end
    end
    
    % Velocity Y
    
    for i=1:NIvt
        for j=1:NJvt
            if (i<=NIue && j>NJue) continue; end
            
            % Top boundaries
            if ( j==NJvt ) V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) ; continue ; end
            if ( (j==NJve) && (i<=NIve) ) V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) ; continue ; end
            
            % Bot boundary
            if ( j==1 ) V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) ; continue ; end
            
            % Right boundary
            if ( i==NIvt ) V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) ; continue ; end
            
            % Left boundary
            if ( i==1 ) V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) ; continue ; end
            
            V ( Vmap (i,j) ) = Vs ( Vmap (i,j) ) - c1 * ( P ( Pmap(i,j) ) - P ( Pmap(i,j-1) ) );
        end
    end
    
    % Saving results        
    if ( mod(nstep-1,nSkip) == 0 )
        nSave = nSave+1;
        PT (:,nSave)= P;
        UT (:,nSave)= U;
        VT (:,nSave)= V;
        
    end
    
end


%% Post-Process



% Plot pressure
figure()

scatter3 (Xp(:,1),Xp(:,2),P,'r','.')
title('pressure');

% Plot velocity x
figure()
%for i=1:1000
scatter3 (Xu(:,1),Xu(:,2),U,'k','.')
title('velocity u');

%end
% Plot velocity y
% figure()
%
% scatter3 (Xv(:,1),Xv(:,2),VT(:,750),'g','.')
% title('velocity v');






























